home *** CD-ROM | disk | FTP | other *** search
/ Turnbull China Bikeride / Turnbull China Bikeride - Disc 2.iso / BARNET / COMPILER / SATHER / !Sather / Library / Math / sa / mat < prev    next >
Text File  |  1996-07-18  |  40KB  |  1,290 lines
  1. ---------------------------> Sather 1.1 source file <--------------------------
  2. -- Copyright (C) International Computer Science Institute, 1994.  COPYRIGHT  --
  3. -- NOTICE: This code is provided "AS IS" WITHOUT ANY WARRANTY and is subject --
  4. -- to the terms of the SATHER LIBRARY GENERAL PUBLIC LICENSE contained in    --
  5. -- the file "Doc/License" of the Sather distribution.  The license is also   --
  6. -- available from ICSI, 1947 Center St., Suite 600, Berkeley CA 94704, USA.  --
  7. --------> Please email comments to "sather-bugs@icsi.berkeley.edu". <----------
  8. -- Author: Kennel <mbk@icsi.berkeley.edu>
  9. -------------------------------------------------------------------
  10. class MAT{ET<$NFE{ET},VT<$VEC{ET,VT}} < $MAT{ET,VT,MAT{ET,VT}} is
  11. -- This is a purely Sather implementation of the generic dense
  12. -- concrete matrix class, with *COLUMN-MAJOR* (fortran-style) 
  13. -- layout.  First index changes most rapidly in stepping thru storage.
  14.  
  15. -- This was written before iterator optimizations - a more natural
  16. -- style using the ind! and other iterators which are now built-in
  17. -- could be adopted. At the time it was written, the compiler
  18. -- did not dod these optimizations.
  19. -- Original note from Matt:n
  20. -- Because iters are not yet optimized or inlined in the compiler I
  21. -- will become a compromise for speed and write out operations with
  22. -- boring conventional loops.  Implementation will tend to use extra
  23. -- index variabless to avoid integer multiplication inside the loop.
  24. -- This may rely on assumed column major ordering.  Also I will
  25. -- attempt to use local variables and not attributes in loops. 
  26.  
  27.  
  28.    private include AREF{ET} aclear->aclear,aget->aget,aset->aset;
  29.  
  30.    -- makes those named features public, others stay private.
  31.  
  32.    readonly attr nr:INT;  -- number of rows
  33.    readonly attr nc:INT;  -- number of columns
  34.    -- built in feature inherited from AREF{*}, asize = nr*nc.
  35.  
  36.    size: INT is return asize end;
  37.    size1:INT is return nr end; 
  38.    size2:INT is return nc end;
  39.    
  40.    element_zero:ET is return ET::zero; end;
  41.    element_one:ET is return ET::one; end; 
  42.    -- These may need be redefined for MATCPX et al.  This is a hack.  Zero
  43.    -- and one should really be features in in ET < $RING_ELT, where
  44.    -- $RING_ELT gives additive and multiplicative identity.
  45.  
  46.    is_same_shape(arg:SAME):BOOL is
  47.       -- useful in preconditions.  Is arg dimensioned the same as self?
  48.       if void(arg) or void(self) then
  49.      return false;
  50.       elsif arg.asize /= asize or arg.nr /= nr then
  51.      return false;
  52.       else return true; end;
  53.    end;
  54.  
  55.    is_same_shape_trans(arg:SAME):BOOL is
  56.       -- does arg^T conform to 'self'?
  57.       if void(arg) or void(self) then
  58.      return false;
  59.       elsif arg.asize /= asize or arg.nr /= nc then
  60.      return false;
  61.       else return true; end;
  62.    end;
  63.  
  64.    fits(arg:SAME):BOOL is
  65.       -- will the contents of 'arg' fit into self, ignoring tranposition.
  66.       if void(arg) or void(self) or arg.asize /= asize then 
  67.      return false; else return true end;
  68.    end;
  69.  
  70.    is_eq(m: SAME): BOOL pre ~void(self) and ~void(m) and is_same_shape(m) is
  71.       loop if elt!/=m.elt! then return false end; end; return true end;
  72.    
  73.    reshape(r,c:INT) pre r*c = asize is
  74.       -- Reshape 'self' to have 'r' rows and 'c' columns but keeping
  75.       -- actual data, as laid out in one dimension, in the same place.
  76.       -- This only changes the bounds, the data will logically move 
  77.       -- rows and columns keeping column major layout.
  78.       -- Example:       [a d]
  79.       --             m1=[b e]
  80.       --                [c f].
  81.       -- m1.nr = 3, m1.nc = 2.   After m1.reshape(2,3):
  82.       --                [a c e]
  83.       --             m1=[b d f]
  84.       --                
  85.       -- This feature is NOT in the abstract supertype because it is
  86.       -- representation dependent.  
  87.       nc := c; nr := r;
  88.    end;
  89.    
  90.    create(r,c:INT):SAME pre (r>0) and (c>0) is
  91.       -- Create a matrix with r rows and c columns
  92.       res::=new(r*c);
  93.       res.nr := r; res.nc := c;
  94.       return(res);
  95.    end;
  96.  
  97.    create(arg:SAME):SAME pre ~void(arg) is
  98.       -- Creates a new matrix with the same dimensions (but
  99.       -- not the same values) as arg
  100.       res::=new(arg.asize);
  101.       res.nr := arg.nr;  res.nc := arg.nc;
  102.       return(res);
  103.    end;
  104.  
  105.    create(a: ARRAY{ARRAY{ET}}): SAME 
  106.    -- Create a new array with the same dimensions and values as
  107.    -- a, which is an array of arrays(rows). 
  108.    -- Assume that all the rows of "a" have the same number of elements
  109.       pre a.size > 0 and a[0].size > 0 is
  110.       sz1 ::= a.size;
  111.       sz2 ::= a[0].size;
  112.       res ::= #SAME(sz1,sz2);
  113.       loop r::=sz1.times!; 
  114.      loop c::=sz2.times!; res[r,c] := a[r][c]; end end;
  115.       return(res);
  116.    end;
  117.    
  118.    str:STR is
  119.       -- The string form of self represented as a list of rows,
  120.       -- eg. "||1.00,2.33|,|4.5,2.8|,|9.7,3.2||".
  121.       sc::=#FSTR + "|";
  122.       loop r::=nr.times!; 
  123.      if r/=0 then sc := sc+ "," end; 
  124.      sc := sc + "|";
  125.      loop c::=nc.times!; 
  126.         if c/=0 then sc := sc+"," end; 
  127.         sc := sc + [r,c].str; end;
  128. --        #OUT+r+","+c+":"+sc+"\n"; end;
  129.      sc := sc+ "|" end; 
  130.       sc := sc+"|"; 
  131.       return(sc.str) end;
  132.  
  133.    dimension_str:STR is
  134.       -- useful for debugging, prints out "nr = self.nr, nc = self.nc"
  135.       return "nr = " + self.nr + ", nc = " + self.nc;
  136.    end;
  137.  
  138.    copy:SAME is
  139.       -- make a value copy.
  140.       res ::= create(self);
  141.       res.inplace_contents(self);
  142.       return(res);
  143.    end;
  144.  
  145.    contents(arg:SAME) is inplace_contents(arg) end;
  146.  
  147.    inplace_contents(arg:SAME) pre is_same_shape(arg) is
  148.       -- copy the contents of arg into self.
  149.       sz ::= asize; i::=0; 
  150.       loop while!(i<sz); self[i] := arg[i];
  151.      i := i+1;
  152.       end;
  153.    end;
  154.  
  155.    inplace_contents_from_function(function:ROUT{INT,INT}:ET) is
  156.       -- Set the contents of the matrix from the function "function"
  157.       loop
  158.      i1 ::= col_ind!;
  159.      loop
  160.         i2 ::= row_ind!;
  161.         [i1,i2] := function.call(i1,i2);
  162.      end;
  163.       end;
  164.    end;
  165.    
  166.    inplace_portion_of_arg(a: SAME) is
  167.       -- Copy into self as much of arg as will fit and return it. Don't
  168.       -- alter other elements.
  169.       nrows: INT := nr.min(a.nr);
  170.       loop r::= nrows.times!;
  171.      loop set_row!(r,a.row_elt!(r)) end 
  172.       end 
  173.    end;
  174.    
  175.    ident:SAME is
  176.       -- Create an identity matrix of the same shape as self
  177.       res ::= create(self);
  178.       -- here we can assume that entries are already set to 0 via create
  179.       d ::= nc.min(nr); 
  180.       i ::= 0; loop while!(i<d);
  181.      res[i,i] := element_one; 
  182.      i := i+1;
  183.       end;
  184.       return res;
  185.    end;
  186.  
  187.    inplace_ident is
  188.       -- Should work for non rectangular matrices too, non
  189.       -- square portion will be set to zero.
  190.       sz ::= asize; i::=0; r::=0; c::=0; lnr ::= nr;
  191.       loop while!(i<sz);
  192.      if (r = c) then [i] := element_one; else [i] := element_zero; end;
  193.      r := r+1;
  194.      if (r >= lnr) then r := 0; c := c+1; end;
  195.      i := i+1;
  196.       end;
  197.    end;
  198.  
  199.    trans:SAME is
  200.       -- Create a new matrix that is the transpose of self
  201.       lnc ::= nc;  res ::= #SAME(lnc,nr);
  202.       -- reverse columns and rows return new value.
  203.       sz ::= asize; j::=0; i::=0;
  204.       loop while!(i<sz);
  205.      res[j] := self[i];
  206.      j := j+lnc;
  207.      diff ::= j-sz; if (diff>=0) then j:=diff+1; end;
  208.      i := i+1;
  209.       end;
  210.       return res;
  211.    end;
  212.  
  213.    inplace_trans is
  214.       -- make self = transpose of self.
  215.       if (nr = nc) then
  216.      -- can do it in place!
  217.      d ::= nr; sz ::= asize;
  218.      i ::= 1; j::=d; c::=0;
  219.      -- Count through lower diagonal indexes and corresponding
  220.      -- upper diagonal indexes in this non-obvious but efficient
  221.      -- manner
  222.      k ::= 0;  -- diagonal index.
  223.      loop
  224.         while!(i < sz);
  225.         tmp ::= [j]; [j] := [i]; [i] := tmp;
  226.         -- flip i and i transpose 
  227.         i := i+1; j:= j+d;
  228.         if j>=sz then
  229.            k := k+d+1;
  230.            i := k+1; -- one down
  231.            j := k+d; -- one to the right
  232.         end;
  233.      end;
  234.       else
  235.      tmp ::= self.trans; -- create temporary.
  236.      self.nc := tmp.nc;
  237.      self.nr := tmp.nr;
  238.      self.inplace_contents(tmp);
  239.       end;
  240.    end;
  241.  
  242.    inplace_arg_trans(arg:SAME) pre fits(arg) is
  243.       -- self <- arg^T
  244.       lnc ::= arg.nc; 
  245.       nc := arg.nr;
  246.       nr := lnc;
  247.       -- reverse columns and rows return new value.
  248.       sz ::= asize; j::=0; i::=0;
  249.       loop while!(i<sz);
  250.      self[j] := arg[i];
  251.      j := j+lnc;
  252.      diff ::= j-sz; if (diff>=0) then j:=diff+1; end;
  253.      i := i+1;
  254.       end;
  255.    end;
  256.  
  257.    times_elt(s:ET):SAME is
  258.       -- Return a new matrix that is self * scalar s
  259.       -- Self is unchanged
  260.       return(scaled_by(s)); 
  261.    end;
  262.  
  263.    inplace(s:ET) is
  264.       -- Set all elements to have the value "s"
  265.       sz ::= asize; i::=0; 
  266.       loop while!(i<sz); [i] := s; i := i+1; end;
  267.    end;
  268.    
  269.    inplace_zero is
  270.       -- Set all elements to zero
  271.       inplace(ET::zero);
  272.    end;
  273.    
  274.    scaled_by(s:ET):SAME is
  275.       res ::= #SAME(self);
  276.       sz ::= asize; i::=0; 
  277.       loop while!(i<sz); res[i] := s*self[i]; i := i+1; end;
  278.       return res;
  279.    end;
  280.  
  281.    inplace_scaled_by(s:ET) is
  282.       sz ::= asize; i::=0; 
  283.       loop while!(i<sz); [i] := s*[i]; i := i+1; end;
  284.    end;
  285.  
  286.    inplace_elements(s:ET) is
  287.       sz ::= asize; i::=0; 
  288.       loop while!(i<sz); [i] := s; i := i+1; end;
  289.    end;
  290.  
  291.    aget(i1,i2:INT):ET pre (i1 >=0) and (i1 <= nr) and (i2>=0) and (i2<=nc) is
  292.       -- The element with indices `[i1,i2]'.
  293.       return([i1+i2*nr]) end;
  294.  
  295.    aset(i1,i2:INT,val:ET) pre (i1 >=0) and (i1 <= nr) and
  296.         (i2>=0) and (i2<=nc) is
  297.       -- Set the element with indices `[i1,i2]' to val.
  298.       [i1+i2*nr]:=val end;      
  299.  
  300.    plus(arg:SAME):SAME is return(plus_arg(arg)); end;
  301.  
  302.    plus_arg(arg:SAME):SAME pre is_same_shape(arg) is
  303.       res ::= #SAME(self);
  304.       res.inplace_arg_plus_arg(self,arg);
  305.       return(res);
  306.    end;
  307.  
  308.    minus(arg:SAME):SAME is return(minus_arg(arg)); end;
  309.  
  310.    minus_arg(arg:SAME):SAME pre is_same_shape(arg) is
  311.       res ::= #SAME(self);
  312.       res.inplace_arg_minus_arg(self,arg);
  313.       return(res);
  314.    end;
  315.  
  316.    inplace_plus_arg(arg:SAME) pre is_same_shape(arg) is
  317.       self.inplace_arg_plus_arg(self,arg);
  318.    end;
  319.  
  320.    inplace_minus_arg(arg:SAME) pre is_same_shape(arg) is
  321.       self.inplace_arg_minus_arg(self,arg);
  322.    end;
  323.  
  324.    plus_arg_trans(arg:SAME):SAME pre is_same_shape_trans(arg) is
  325.       res ::= #SAME(self);
  326.       res.inplace_arg_plus_arg_trans(self,arg);
  327.       return res;
  328.    end;
  329.  
  330.    minus_arg_trans(arg:SAME):SAME pre is_same_shape_trans(arg) is
  331.       res ::= #SAME(self);
  332.       res.inplace_arg_minus_arg_trans(self,arg);
  333.       return res;
  334.    end;
  335.  
  336.    inplace_plus_arg_trans(arg:SAME) pre is_same_shape_trans(arg) is
  337.       self.inplace_arg_plus_arg_trans(self,arg);
  338.    end;
  339.  
  340.    inplace_minus_arg_trans(arg:SAME) pre is_same_shape_trans(arg) is
  341.       self.inplace_arg_minus_arg_trans(self,arg);
  342.    end;
  343.  
  344. -- Three way operations: self <- arg1^{T} (+/-) arg2{^T}
  345. -- can write others in terms of these.
  346.  
  347.    inplace_arg_plus_arg(arg1,arg2:SAME) pre is_same_shape(arg1) and 
  348.         arg1.is_same_shape(arg2) is
  349.       sz ::= asize; i::=0; 
  350.       loop while!(i<sz); self[i] := arg1[i] + arg2[i];
  351.      i := i+1;
  352.       end;
  353.    end;
  354.  
  355.    inplace_arg_minus_arg(arg1,arg2:SAME) pre is_same_shape(arg1) and 
  356.         arg1.is_same_shape(arg2) is
  357.       sz ::= asize; i::=0; 
  358.       loop while!(i<sz); self[i] := arg1[i] - arg2[i];
  359.      i := i+1;
  360.       end;
  361.    end;
  362.  
  363.    inplace_arg_plus_arg_trans(arg1,arg2:SAME) pre fits(arg1)
  364.         and arg1.is_same_shape_trans(arg2) and ~SYS::ob_eq(self,arg2) is
  365.       lnc ::= arg1.nc;
  366.       nc := lnc;  nr := arg1.nr;
  367.       sz ::= asize; j::=0; i::=0;
  368.       loop while!(i<sz);
  369.      self[i] := arg1[i] + arg2[j];
  370.      j := j+lnc;
  371.      diff ::= j-sz; if (diff>=0) then j:=diff+1; end;
  372.      i := i+1;
  373.       end;
  374.    end;
  375.  
  376.    inplace_arg_minus_arg_trans(arg1,arg2:SAME) pre fits(arg1)
  377.         and arg1.is_same_shape_trans(arg2) and ~SYS::ob_eq(self,arg2) is
  378.       lnc ::= arg1.nc;
  379.       nc := lnc;  nr := arg1.nr;
  380.       sz ::= asize; j::=0; i::=0;
  381.       loop while!(i<sz);
  382.      self[i] := arg1[i] - arg2[j];
  383.      j := j+lnc;
  384.      diff ::= j-sz; if (diff>=0) then j:=diff+1; end;
  385.      i := i+1;
  386.       end;
  387.    end;
  388.  
  389.    inplace_arg_trans_plus_arg_trans(arg1,arg2:SAME) pre fits(arg1) 
  390.      and arg1.is_same_shape(arg2) and
  391.         ~SYS::ob_eq(self,arg1) and ~SYS::ob_eq(self,arg2) is
  392.       lnc ::= arg1.nc;
  393.       nr := lnc;  nc := arg1.nr;
  394.       sz ::= asize; j::=0; i::=0;
  395.       loop while!(i<sz);
  396.      self[j] := arg1[i] + arg2[i];
  397.      j := j+lnc;
  398.      diff ::= j-sz; if (diff>=0) then j:=diff+1; end;
  399.      i := i+1;
  400.       end;
  401.    end;
  402.  
  403.    inplace_arg_trans_minus_arg_trans(arg1,arg2:SAME) pre fits(arg1)
  404.      and arg1.is_same_shape(arg2) and
  405.         ~SYS::ob_eq(self,arg1) and ~SYS::ob_eq(self,arg2) is
  406.       lnc ::= arg1.nc;
  407.       nr := lnc;  nc := arg1.nr;
  408.       sz ::= asize; j::=0; i::=0;
  409.       loop while!(i<sz);
  410.      self[j] := arg1[i] - arg2[i];
  411.      j := j+lnc;
  412.      diff ::= j-sz; if (diff>=0) then j:=diff+1; end;
  413.      i := i+1;
  414.       end;
  415.    end;
  416.  
  417. -- with scaling, euphemisms for elementary 3 way operations.
  418.    plus_scaled_arg(s:ET,arg:SAME):SAME is
  419.       res ::= #SAME(self);
  420.       res.inplace_arg_plus_scaled_arg(self,s,arg);
  421.       return res;
  422.    end;
  423.  
  424.    inplace_plus_scaled_arg(s:ET,arg:SAME) is
  425.       self.inplace_arg_plus_scaled_arg(self,s,arg);
  426.    end;
  427.  
  428.    plus_scaled_arg_trans(s:ET,arg:SAME):SAME is
  429.       res ::= #SAME(self);
  430.       res.inplace_arg_plus_scaled_arg_trans(self,s,arg);
  431.       return(res);
  432.    end;
  433.  
  434.    inplace_plus_scaled_arg_trans(s:ET,arg:SAME) is
  435.       self.inplace_arg_plus_scaled_arg_trans(self,s,arg);
  436.    end;
  437.  
  438. -- with scaling, 3 way operations.
  439.    inplace_arg_plus_scaled_arg(arg1:SAME,s:ET,arg2:SAME) pre
  440.         is_same_shape(arg1) and arg1.is_same_shape(arg2) is
  441.       sz ::= asize; i::=0; 
  442.       loop while!(i<sz); self[i] := arg1[i] + s*arg2[i];
  443.      i := i+1;
  444.       end;
  445.    end;
  446.  
  447.    inplace_arg_plus_scaled_arg_trans(arg1:SAME,s:ET,arg2:SAME) pre fits(arg1)
  448.         and arg1.is_same_shape_trans(arg2) and ~SYS::ob_eq(self,arg2) is
  449.       lnc ::= arg1.nc;
  450.       nc := lnc;  nr := arg1.nr;
  451.       sz ::= asize; j::=0; i::=0;
  452.       loop while!(i<sz);
  453.      self[i] := arg1[i] + s*arg2[j];
  454.      j := j+lnc;
  455.      diff ::= j-sz; if (diff>=0) then j:=diff+1; end;
  456.      i := i+1;
  457.       end;
  458.    end;
  459.    
  460. --
  461. -- Matrix Multiplication.
  462. --
  463.    times(arg:SAME):SAME pre nc = arg.nr is
  464.       res ::= #SAME(nr,arg.nc);
  465.       res.inplace_arg_times_arg(self,arg);
  466.       return res;
  467.    end;
  468.  
  469.    trans_times_arg(arg:SAME):SAME pre nr = arg.nr is
  470.       res ::= #SAME(nc,arg.nc);
  471.       res.inplace_arg_trans_times_arg(self,arg);
  472.       return res;
  473.    end;
  474.  
  475.    times_arg_trans(arg:SAME):SAME pre nc = arg.nc is
  476.       res ::= #SAME(nr,arg.nr);
  477.       res.inplace_arg_times_arg_trans(self,arg);
  478.       return res;
  479.    end;
  480.  
  481.    trans_times_arg_trans(arg:SAME):SAME pre nr = arg.nc is
  482.       res ::= #SAME(nc,arg.nr);
  483.       res.inplace_arg_trans_times_arg_trans(self,arg);
  484.       return res;
  485.    end;
  486.  
  487. -- in place versions:
  488.  
  489.    inplace_arg_times_arg(arg1,arg2:SAME) pre nr = arg1.nr and
  490.         nc = arg2.nc and arg1.nc = arg2.nr is
  491.       -- self := arg1 * arg2.
  492.       -- For all i,j, self[i,j] = sum(k) arg1[i,k] * arg2[k,j]
  493.  
  494.       selfnr ::= nr; selfnc ::= nc;
  495.       a1nr ::= arg1.nr; a1nc ::= arg1.nc;
  496.       a2nr ::= arg2.nr; a2nc ::= arg2.nc;
  497.       -- This is what you need to do to get reasonable performance on
  498.       -- matrix codes: take out of the loop the array index multiplication,
  499.       -- i.e. the definition of [i,j] ->  [i+j*nr]
  500.       -- Eventually the compiler will figure this out.  Note that
  501.       -- this is not as trivial as you might disdainfully think, as in
  502.       -- Sather, and unlike Fortran, the size of the array could conceivably
  503.       -- change during a loop through any call that might modify nr or nc.
  504.       -- The array reference itself might change as well.   Sather 1.0.8 and
  505.       -- above ought to have dataflow optimizations that can do this in
  506.       -- most cases, but I'd rather be safe than sorry here.
  507.       -- 
  508.       -- The following routine seems to generate assembly code identical
  509.       -- to a naive C language implementation for its inner loops with GNU
  510.       -- CC for ET=FLT. 
  511.  
  512.       i::=0;
  513.       loop while!(i<selfnr);
  514.      j::=0;
  515.      loop while!(j<selfnc);
  516.         sum ::= element_zero;
  517.         k::=0;
  518.         loop while!(k < a1nc); 
  519.            sum := sum+arg1[i+k*a1nr]*arg2[k+j*a2nr]; k := k+1;
  520.         end;
  521.         [i+j*selfnr] := sum;
  522.         j := j+1;
  523.      end;
  524.      i := i+1;
  525.       end;
  526.    end;
  527.  
  528.    inplace_arg_times_arg_fast(arg1,arg2:SAME) 
  529.    -- self := arg1 * arg2.
  530.    -- Basic algorithm:
  531.    -- For all i,j, self(i,j) = sum(k) arg1[i,k] * arg2[k,j]
  532.    -- t1 = i+j*nr, index through "self".  0 .. asize-1
  533.    -- t2 = i+arg1.nr*k                    0 .. arg1.asize-1
  534.    -- t3 = k+j*arg2.nr
  535.       pre nr = arg1.nr and  nc = arg2.nc and arg1.nc = arg2.nr is
  536.       lnr ::= nr; lnc ::= nc;
  537.       sz ::= asize;  szarg1 ::= arg1.asize;
  538.       lk ::= arg1.nc;  -- lk is the bounds of the innermost loop.
  539.  
  540.       j::=0;
  541.       idx0 ::= 0;
  542.       idx2init ::= 0;
  543.       loop while!(j<lnc);
  544.      i::=0;
  545.      loop while!(i<lnr);
  546.         sum ::= element_zero;
  547.  
  548.         idx1 ::= i;
  549.         idx2 ::= idx2init; -- was "j*lk"
  550.         assert(idx2init = j*lk);
  551.         k::=lk-1; -- do innermost loop lk times.
  552.         loop while!(k >= 0);
  553.            sum := sum+arg1[idx1]*arg2[idx2];
  554.            idx1 := idx1+lnr; -- lnr is = arg1.nr
  555.            idx2 := idx2+1;
  556.            k := k-1;
  557.         end;
  558. --        [i+j*nr] := sum;
  559.         [idx0] := sum;
  560.         assert (i+j*lnr = idx0);
  561.         idx0 := idx0+1;
  562.         i := i+1;
  563.      end;
  564.      j := j+1;
  565.      idx2init := idx2init + lk;
  566.       end;
  567.    end;
  568.  
  569.    inplace_arg_trans_times_arg(arg1,arg2:SAME) 
  570.    -- self := arg1^T * arg2.
  571.    -- For all i,j, self[i,j] = sum(k) arg1[k,i] * arg2[k,j]
  572.    -- Likely faster than arg * arg, as inner loop is unit stride
  573.    -- for both arrays.
  574.       pre nr = arg1.nc and nc = arg2.nc and arg1.nr = arg2.nr 
  575.    is
  576.       selfnr ::= nr; selfnc ::= nc;
  577.       a1nr ::= arg1.nr; a1nc ::= arg1.nc;
  578.       a2nr ::= arg2.nr; a2nc ::= arg2.nc;
  579.  
  580.       i::=0;
  581.       loop while!(i<selfnr);
  582.      j::=0;
  583.      loop while!(j<selfnc);
  584.         sum ::= element_zero;
  585.         k::=0;
  586.         loop while!(k < a1nr); 
  587.            sum := sum+arg1[k+i*a1nr]*arg2[k+j*a2nr]; k := k+1;
  588.         end;
  589.         [i+j*selfnr] := sum;
  590.         j := j+1;
  591.      end;
  592.      i := i+1;
  593.       end;
  594.    end;
  595.    
  596.    inplace_arg_times_arg_trans(arg1,arg2:SAME) 
  597.    -- self := arg1 * arg2^T.
  598.    -- For all i,j, self[i,j] = sum(k) arg1[i,k] * arg2[j,k]
  599.       pre nr = arg1.nr and nc = arg2.nr and arg1.nc = arg2.nc 
  600.    is
  601.       selfnr ::= nr; selfnc ::= nc;
  602.       a1nr ::= arg1.nr; a1nc ::= arg1.nc;
  603.       a2nr ::= arg2.nr; a2nc ::= arg2.nc;
  604.  
  605.       i::=0;
  606.       loop while!(i<selfnr);
  607.      j::=0;
  608.      loop while!(j<selfnc);
  609.         sum ::= element_zero;
  610.         k::=0;
  611.         loop while!(k < a1nc);  -- would be better to rearrange loop
  612.            sum := sum+arg1[i+k*a1nr]*arg2[j+k*a2nr]; k := k+1;
  613.         end;
  614.         [i+j*selfnr] := sum;
  615.         j := j+1;
  616.      end;
  617.      i := i+1;
  618.       end;
  619.    end;
  620.  
  621.    inplace_arg_trans_times_arg_trans(arg1,arg2:SAME) 
  622.    -- self := arg1^T * arg2^T.
  623.    -- For all i,j, self[i,j] = sum(k) arg1[k,i] * arg2[j,k]
  624.       pre nr = arg1.nc and nc = arg2.nr and arg1.nr = arg2.nc 
  625.    is
  626.       selfnr ::= nr; selfnc ::= nc;
  627.       a1nr ::= arg1.nr; a1nc ::= arg1.nc;
  628.       a2nr ::= arg2.nr; a2nc ::= arg2.nc;
  629.  
  630.       i::=0;
  631.       loop while!(i<selfnr);
  632.      j::=0;
  633.      loop while!(j<selfnc);
  634.         sum ::= element_zero;
  635.         k::=0;
  636.         loop while!(k < a1nr); 
  637.            sum := sum+arg1[k+i*a1nr]*arg2[j+k*a2nr]; k := k+1;
  638.         end;
  639.         [i+j*selfnr] := sum;
  640.         j := j+1;
  641.      end;
  642.      i := i+1;
  643.       end;
  644.    end;
  645.    
  646.    inplace_times_diagonal(v1:VT) is
  647.       -- Set self to be the product of itself with a diagonal matrix whose
  648.       -- diagonal entries are the components of v1 truncated or extended
  649.       -- with zeroes to be the correct size.
  650.       -- Scale the columns of self with the elements of v1
  651.       -- Self <- Self*v1 (as diagonal of a matrix)
  652.       loop c::=nc.min(v1.dim).times!;
  653.      loop set_col!(c,v1[c]*col_elt!(c)) end 
  654.       end;
  655.       loop c::=v1.dim.upto!(nc-1); 
  656.      loop set_col!(c,ET::zero)  end; 
  657.       end;
  658.    end;
  659. --
  660. -- MATRIX/VECTOR operations
  661. --
  662. -- create a new argument
  663.  
  664.    times_vec(arg:VT):VT pre arg.dim = self.nc is
  665.       -- This is a syntactical sugar expressions (m*v) so it
  666.       -- doesn't follow the naming convention.
  667.       res:VT := #VT(nr);
  668.       times_vec_into_vec(arg,res);
  669.       return res;
  670.    end;
  671.    
  672.    trans_times_vec(arg:VT):VT pre arg.dim = self.nr is
  673.       res:VT := #VT(nc);
  674.       trans_times_vec_into_vec(arg,res);
  675.       return res;
  676.    end;
  677.    
  678. -- in place
  679.  
  680.    times_vec_into_vec(arg,dest:VT) 
  681.       pre arg.dim = self.nc and dest.dim = self.nr 
  682.    is
  683.       selfnr ::= nr; selfnc ::= nc;
  684.       i ::= 0;
  685.       loop while!(i<selfnr);
  686.      sum ::= element_zero;
  687.      j::= 0;
  688.      loop while!(j<selfnc);
  689.         sum := sum + self[i,j]*arg[j];
  690.         j := j+1;
  691.      end;
  692.      dest[i] := sum;
  693.      i:=i+1;
  694.       end;
  695.    end;
  696.  
  697.    trans_times_vec_into_vec(arg,dest:VT) 
  698.       pre arg.dim = self.nr and dest.dim = self.nc 
  699.    is
  700.       selfnr ::= nr; selfnc ::= nc;
  701.       i ::= 0;
  702.       loop while!(i<selfnc);
  703.      sum ::= element_zero;
  704.      j::= 0;
  705.      loop while!(j<selfnr);
  706.         sum := sum + self[j,i]*arg[j];
  707.         j := j+1;
  708.      end;
  709.      dest[i] := sum;
  710.      i:=i+1;
  711.       end;
  712.    end;
  713.  
  714.    times_scaled_vec_into_vec(s:ET,arg,dest:VT) 
  715.       pre arg.dim = self.nc and dest.dim = self.nr 
  716.    is
  717.       selfnr ::= nr; selfnc ::= nc;
  718.       i ::= 0;
  719.       loop while!(i<selfnr);
  720.      sum ::= element_zero;
  721.      j::= 0;
  722.      loop while!(j<selfnc);
  723.         sum := sum + self[i,j]*arg[j];
  724.         j := j+1;
  725.      end;
  726.      dest[i] := s*sum;
  727.      i:=i+1;
  728.       end;
  729.    end;
  730.  
  731.    trans_times_scaled_vec_into_vec(s:ET,arg,dest:VT) 
  732.       pre arg.dim = self.nr and dest.dim = self.nc 
  733.    is
  734.       selfnr ::= nr; selfnc ::= nc;
  735.       i ::= 0;
  736.       loop while!(i<selfnc);
  737.      sum ::= element_zero;
  738.      j::= 0;
  739.      loop while!(j<selfnr);
  740.         sum := sum + self[j,i]*arg[j];
  741.         j := j+1;
  742.      end;
  743.      dest[i] := s*sum;
  744.      i:=i+1;
  745.       end;
  746.    end;
  747.  
  748.    inplace_plus_scaled_vec_times_vec(s:ET,v1,v2:VT) 
  749.       -- self := self + s*v1*v2^T
  750.       -- (Add scaled outer product of v1 and v2 to self.  A BLAS operation)
  751.       --
  752.       pre v1.dim = nr and v2.dim = nc 
  753.    is
  754.          selfnr ::= nr; selfnc ::= nc;
  755.       i ::= 0;
  756.       loop while!(i<selfnc);
  757.      j::= 0;
  758.      loop while!(j<selfnr);
  759.         [j,i] := [j,i] +s*v1[j]*v2[i];
  760.         j := j+1;
  761.      end;
  762.      i:=i+1;
  763.       end;
  764.    end;
  765.    
  766. --
  767. -- MATRIX/VECTOR MANIPULATION
  768. --
  769.    create_col_matrix(arg:VT):SAME is
  770.       d ::= arg.dim;
  771.       res ::= create(d,1);
  772.       i ::= 0; loop while!(i<d); res[i] := arg[i]; i:=i+1; end; 
  773.    end;
  774.  
  775.    create_row_matrix(arg:VT):SAME is
  776.       d ::= arg.dim;
  777.       res ::= create(1,d);
  778.       i ::= 0; loop while!(i<d); res[i] := arg[i]; i := i+1; end; 
  779.    end;
  780.  
  781.    col(i:INT):VT is
  782.       selfnr ::= nr;
  783.       res:VT := #VT(selfnr);
  784.       j::=0; loop while!(j<selfnr); res[j] := [j,i]; j := j+1; end;
  785.    end;
  786.  
  787.    row(i:INT):VT is
  788.       selfnc ::= nc;
  789.       res:VT := #VT(selfnc);
  790.       j::=0; loop while!(j<selfnc); res[j] := [i,j]; j := j+1; end;
  791.    end;
  792.  
  793.    col(i:INT,v:VT) pre v.dim = nr is
  794.       d ::= v.dim;
  795.       j ::= 0;
  796.       loop while!(j < d);  [j,i] := v[j]; j:= j+1; end;
  797.    end;
  798.  
  799.    row(i:INT,v:VT) pre v.dim = nc is
  800.       d ::= v.dim;
  801.       j ::= 0;
  802.       loop while!(j < d);  [i,j] := v[j]; j:= j+1; end;
  803.    end;
  804.  
  805.    inplace_scaled_col(s:ET,i:INT) pre i.is_bet(0,nc-1) is 
  806.       d ::= nr; j ::= 0;
  807.       loop while!(j < d);  [j,i] := [j,i]*s; j:= j+1; end;
  808.    end;
  809.    
  810.    inplace_scaled_row(s:ET,i:INT) pre i.is_bet(0,nr-1) is
  811.       d ::= nc; j ::= 0;
  812.       loop while!(j < d);  [i,j] := [i,j]*s; j:= j+1; end;
  813.    end;
  814.    
  815.    inplace_col_plus_scaled_vec(i:INT,s:ET,v:VT) 
  816.       pre i.is_bet(0,nc-1) and v.dim = nr 
  817.    is
  818.       d ::= nr; j ::= 0;
  819.       loop while!(j < d);  [j,i] := [j,i] + s*v[j]; j:= j+1; end;
  820.    end;
  821.  
  822.    inplace_row(i: INT,v: VT) pre i.is_bet(0,nr-1) and v.dim = nc 
  823.    is
  824.       d ::= nc; j ::= 0;
  825.       loop while!(j < d);  [i,j] := v[j]; j:= j+1; end;
  826.    end;
  827.    
  828.    inplace_row_plus_scaled_vec(i:INT,s:ET,v:VT) 
  829.       pre i.is_bet(0,nr-1) and v.dim = nc is
  830.       d ::= nc; j ::= 0;
  831.       loop while!(j < d);  [i,j] := [i,j] + s*v[j]; j:= j+1; end;
  832.    end;
  833.    
  834.    inplace_swapped_col(i:INT,v:VT) pre i.is_bet(0,nc-1) and v.dim = nr is
  835.       d ::= nr; j ::= 0;
  836.       loop while!(j < d);  t ::= [j,i]; [j,i] := v[j]; v[j] := t; j:= j+1; end;
  837.    end;
  838.    
  839.    inplace_swapped_row(i:INT,v:VT) pre i.is_bet(0,nr-1) and v.dim = nc is
  840.       d ::= nc; j ::= 0;
  841.       loop while!(j < d);  t ::= [i,j]; [i,j] := v[j]; v[j] := t; j:= j+1; end;
  842.    end;
  843.    
  844.    submatrix(lr,ur:INT, lc,uc:INT):SAME pre 0 <= lr and lr <= ur and
  845.         ur < nr and 0 <= lc and lc <= uc and uc < nr is
  846.       -- return a submatrix
  847.       rnr ::= ur-lr+1; rnc ::= uc-lc+1;
  848.       res ::= #SAME(rnr,rnc);
  849.       j::=0;
  850.       loop
  851.      while!(j<rnc);
  852.      i ::= 0;
  853.      loop
  854.         while!(i<rnr); res[i,j] := [i+lr,j+lc]; 
  855.         i := i+1;
  856.      end;
  857.      j := j+1;
  858.       end;
  859.       return res;
  860.    end; 
  861.      
  862.    inplace_submatrix_to_arg(lr,ur:INT, lc,uc:INT,arg:SAME) 
  863.    -- Set the submatrix of self given by [lr..ur,lc..uc] to be arg.
  864.    -- with proposed syntax extension:
  865.    --   m.submatrix(1,3,2,4) := m2;
  866.       pre 0 <= lr and lr <= ur and ur < nr and 0 <= lc 
  867.      and lc <= uc and uc < nr and arg.nr = ur-nr+1 and arg.nc = uc-lc+1 
  868.    is
  869.       anr ::= arg.nr; anc ::= arg.nc;
  870.       j::=0;
  871.       loop
  872.      while!(j<anc);
  873.      i ::= 0;
  874.      loop
  875.         while!(i<anr); [i+lr,j+lc] := arg[i,j];
  876.         i := i+1;
  877.      end;
  878.      j := j+1;
  879.       end;
  880.    end;     
  881.  
  882.    inplace_swapped(arg:SAME) pre is_same_shape(arg) is 
  883.       -- Swap the elements of "arg" with self
  884.       sz ::= asize;
  885.       i ::= 0;loop while!(i<sz); t ::= arg[i]; arg[i] := [i]; [i] := t; i:=i+1; end;
  886.    end; 
  887.  
  888.    ind1!: INT is
  889.       -- Yield each value of the first index in order. The rows
  890.       loop yield(size1.times!); end 
  891.    end;
  892.    
  893.    ind2!:INT is
  894.       -- Yield each value of the second index in order. The columns
  895.       loop yield(size2.times!); end 
  896.    end;
  897.  
  898.    row_ind!: INT is
  899.       -- Yield each value of the first index in order. The rows
  900.       loop yield(size1.times!); end 
  901.    end;
  902.    
  903.    col_ind!:INT is
  904.       -- Yield each value of the second index in order. The columns
  905.       loop yield(size2.times!); end 
  906.    end;
  907.    
  908.    inds!:TUP{INT,INT} is
  909.       -- Yield tuples of the indices of self in lexicographical order.
  910.       loop row ::=size1.times!;
  911.      loop yield(#TUP{INT,INT}(row,size2.times!)); end 
  912.       end 
  913.    end;
  914.  
  915.    elt!: ET pre ~void(self) is
  916.       -- Yield all elements in row major order
  917.       loop yield(aelt!) end; 
  918.    end;
  919.    
  920.    inplace_elt!(val:ET) is
  921.       -- Set all elements in row major order
  922.       loop aset!(val); yield; end 
  923.    end;
  924.  
  925.  
  926.    row_elt!(once row:INT):ET pre ~void(self) is
  927.       -- Yield elements by varying index 2 and holding index 1 at `row'.
  928.       -- The elements of a row "row"
  929.       loop yield(aelt!(row,size2,size1)); end 
  930.    end;
  931.  
  932.    col_elt!(once col:INT):ET pre ~void(self) is
  933.       -- Yield elements by varying index 1 and holding index 2 at `col'.
  934.       -- The elements of a "column" col
  935.       loop yield(aelt!(col*size1,size1,1)); end 
  936.    end;
  937.  
  938.    set_row!(once i1:INT, val:ET) pre ~void(self) is
  939.       -- Set to val elements with varying index 2 and index 1 fixed at `i1'.
  940.       -- i.e. setting the row i1
  941.       loop aset!(i1,size2,size1,val); yield end 
  942.    end;
  943.  
  944.    set_col!(once i2:INT, val:ET) pre ~void(self) is
  945.       -- Set to val elements with varying index 1 and index 2 fixed at `i2'.
  946.       -- i.e. setting the column i2 
  947.       loop aset!(i2*size1,size1,val); yield end 
  948.    end;
  949.  
  950.    inplace_row!(once row:INT, val:ET) pre ~void(self) is
  951.       -- Set to val elements with varying index 2 and index 1 fixed at `row'.
  952.       -- i.e. setting a row "row"
  953.       loop aset!(row,size2,size1,val); yield end end;
  954.  
  955.    inplace_col!(once col:INT, val:ET) pre ~void(self) is
  956.       -- Set to val elements with varying index 1 and index 2 fixed at `col'.
  957.       -- i.e. setting the column col 
  958.       loop aset!(col*size1,size1,val); yield end 
  959.    end;
  960.    
  961.    diag_elt!: ET pre ~void(self) is
  962.       -- Yield values along the diagonal (square in smaller dimension)
  963.       loop ind ::= (size1.min(size2)).times!; yield([ind,ind]) end; 
  964.    end;
  965.  
  966.    inplace_diag_elt!(val:ET) pre ~void(self) is
  967.       -- Set values along the diagonal (square in smaller dimension)
  968.       loop id ::= (size1.min(size2)).times!; [id,id] := val; yield;  end; 
  969.    end;
  970.  
  971.    elt1!(once i1:INT):ET pre ~void(self) is
  972.       -- Yield elements by varying index 2 and holding index 1 at `i1'.
  973.       -- The elements of a row "i1"
  974.       -- this is the same as row_elt!
  975.       loop yield(aelt!(i1,size2,size1)); end 
  976.    end;
  977.  
  978.    elt2!(once i2:INT):ET pre ~void(self) is
  979.       -- Yield elements by varying index 1 and holding index 2 at `i2'.
  980.       -- The elements of a "column" i2
  981.       -- this is the same as col_elt!
  982.       loop yield(aelt!(i2*size1,size1,1)); end 
  983.    end;
  984.    
  985.    set1!(once i1:INT, val:ET) pre ~void(self) is
  986.       -- Set to val elements with varying index 2 and index 1 fixed at `i1'.
  987.       -- i.e. setting the row i1
  988.       -- this is the same as set_row!
  989.       loop aset!(i1,size2,size1,val); yield end 
  990.    end;
  991.  
  992.    set2!(once i2:INT, val:ET) pre ~void(self) is
  993.       -- Set to val elements with varying index 1 and index 2 fixed at `i2'.
  994.       -- i.e. setting the column i2 
  995.       -- this is the same as set_col!
  996.       loop aset!(i2*size1,size1,val); yield end 
  997.    end;
  998.  
  999.    
  1000. end; 
  1001. -----------------------------------------------------------------------
  1002. class NUMERIC_MAT{ET<$NUMBER{ET},VT<$VEC{ET,VT}} < $MAT{ET,VT,NUMERIC_MAT{ET,VT}} is
  1003.    -- Functions that don't work on complex numbers but do work on 
  1004.    -- reals
  1005.    include MAT{ET,VT};
  1006.    
  1007.    destructive_invert: SAME pre nr=nc is
  1008.       -- Destructive invert self, return a new matrix.  Similar to the
  1009.       -- Gaussian algorithm.  Raise "Division by Zero" for singular
  1010.       -- matrices.  The Gaussian algorithm is from Sedgewick:
  1011.       -- "Algorithms", pp 57-65.
  1012.  
  1013.       -- eliminate
  1014.       A::=create(nr,nr);
  1015.       A.inplace_ident;
  1016.       loop i::=0.upto!(nr-1);
  1017.          max::=i;
  1018.          loop j::=(i+1).upto!(nr-1);
  1019.             if [j,i].abs > [max,i].abs then max:=j end;
  1020.          end; -- loop
  1021.          loop k::=i.upto!(nr-1);
  1022.             t::=[i,k]; [i,k]:=[max,k]; [max,k]:=t;
  1023.          end; -- loop
  1024.      loop
  1025.         t::=A.row_elt!(i);
  1026.         A.set_row!(i,A.row_elt!(max));
  1027.         A.set_row!(max,t);
  1028.      end; -- loop
  1029.          loop j::=(i+1).upto!(nr-1);
  1030.         loop
  1031.            A.set_row!(j, A.row_elt!(j) - A.row_elt!(i)*[j,i]/[i,i]);
  1032.         end; -- loop
  1033.         loop k::=(nr-1).downto!(i);
  1034.            [j,k] := [j,k] - [i,k]*[j,i]/[i,i];
  1035.         end; -- loop
  1036.          end; -- loop
  1037.       end; -- loop
  1038.       -- substitute
  1039.       loop j::=(nr-1).downto!(0);
  1040.      loop col::=A.col_ind!;
  1041.         t::=ET::zero;
  1042.         loop k::=(j+1).upto!(nr-1);
  1043.            t := t + [j,k]*A[col,k];
  1044.         end; -- loop
  1045.         A[j,col] := (A[j,col]-t)/[j,j];
  1046.      end; -- loop
  1047.       end; -- loop
  1048.       return A
  1049.    end; -- ninv
  1050.    
  1051.    invert: SAME pre nr=nc is return copy.destructive_invert end;
  1052.    -- Non-destructive inversion
  1053.    
  1054.    destructive_det: ET pre nr=nc is
  1055.       -- Destructive determinant of self.  Similar to the Gaussian algorithm.
  1056.       -- The Gaussian algorithm is from Sedgewick: "Algorithms", pp 57-65.
  1057.       
  1058.       -- eliminate
  1059.       loop i::=0.upto!(nr-1);
  1060.          max::=i;
  1061.          loop j::=(i+1).upto!(nr-1);
  1062.             if [j,i].abs > [max,i].abs then max:=j end;
  1063.          end; -- loop
  1064.          loop k::=i.upto!(nr-1);
  1065.             t::=[i,k]; [i,k]:=[max,k]; [max,k]:=t;
  1066.          end; -- loop
  1067.          loop j::=(i+1).upto!(nr-1);
  1068.         loop k::=(nr-1).downto!(i);
  1069.            [j,k] := [j,k] - [i,k]*[j,i]/[i,i];
  1070.         end; -- loop
  1071.          end; -- loop
  1072.       end; -- loop
  1073.       res::=ET::one; loop res:=res*diag_elt!; end;
  1074.       return res
  1075.    end; -- destructive_det
  1076.    
  1077.    det: ET pre nr=nc is return copy.destructive_det end;
  1078.    -- Non destructive determinant routine
  1079.    
  1080.    destructive_gauss(x:VT) pre nc=nr and x.dim=nc is
  1081.       -- Solve linear equations.  Result in x.  self is changed.
  1082.       -- Raise "Division by Zero" for singular matrices.
  1083.       -- This algorithm is from Sedgewick: "Algorithms", pp 57-65.
  1084.       
  1085.       -- eliminate
  1086.       loop i::=0.upto!(nc-1);
  1087.          max::=i;
  1088.          loop j::=(i+1).upto!(nc-1);
  1089.             if [j,i].abs > [max,i].abs then max:=j end;
  1090.          end; -- loop
  1091.          loop k::=i.upto!(nc-1);
  1092.             t::=[i,k]; [i,k]:=[max,k]; [max,k]:=t;
  1093.          end; -- loop
  1094.          t::=x[i]; x[i]:=x[max]; x[max]:=t;
  1095.          loop j::=(i+1).upto!(nc-1);
  1096.             x[j] := x[j] - x[i]*[j,i]/[i,i];
  1097.             loop k::=(nc-1).downto!(i);
  1098.                [j,k] := [j,k] - [i,k]*[j,i]/[i,i];
  1099.             end; -- loop
  1100.          end; -- loop
  1101.       end; -- loop
  1102.       -- substitute
  1103.       loop j::=(nc-1).downto!(0);
  1104.      t::=ET::zero;
  1105.          loop k::=(j+1).upto!(nc-1);
  1106.             t := t + [j,k]*x[k];
  1107.          end; -- loop
  1108.          x[j] := (x[j]-t)/[j,j];
  1109.       end; -- loop
  1110.    end; -- destructive_gauss
  1111.    
  1112.    gauss(x:VT):VT pre nc=nr and x.dim=nc is
  1113.       -- Non destructive gaussian routine
  1114.       y::=x.copy; copy.destructive_gauss(y); return y 
  1115.    end;
  1116.  
  1117. end;
  1118. ----------------------------------------------------------------------------
  1119. class MAT < $MAT{FLT,VEC,MAT} is
  1120.    -- Includes some functions that only work with FLTs
  1121.    -- Generalizing these functions is possible, but would require definitions
  1122.    -- of machine epsilon in the numeric classes
  1123.    include NUMERIC_MAT{FLT,VEC};
  1124.    
  1125.    inplace_uniform_random is 
  1126.       -- Become self's entries uniform in `[0.,1.)' 
  1127.       loop r::=row_ind!; loop [r,col_ind!] := RND::uniform.flt; end; end;
  1128.    end; -- inplace_uniform_random
  1129.  
  1130.    svd_in(a:MAT, w:VEC, v:MAT)
  1131.    -- Computes the singular value decomposition of `self = a w v^T'.
  1132.    -- `a' must be `max(nr,nc)' by `nc', `w' length `nc', `v' is `nc' by
  1133.    -- `nc'. `Self' is unchanged, `a', `w', `v' are altered.
  1134.       pre a.nr=nr.max(nc) and a.nc=nc and 
  1135.      w.dim=nc and v.nr=nc and v.nc=nc  is
  1136.       -- fill in a with self and extra zero rows if necessary
  1137.       a.inplace_zero; a.inplace_portion_of_arg(self); 
  1138.       NR_SVD::svd(a,w,v);    -- Start with Numerical Recipes version.
  1139.      -- Eventually use a better algorithm.
  1140.    end; -- svd_in
  1141.    
  1142.    svd_back_sub(u:MAT, w:VEC, v:MAT, b,x:VEC)
  1143.    -- Solves `a.x=b' for `x' when `a=u.d.v^T' is the svd of `a'.
  1144.       pre u.nc=w.dim and v.nr=u.nc and v.nc=u.nc and
  1145.      b.dim=u.nr and x.dim=u.nc is 
  1146.  
  1147.       tmp::=#VEC(u.nc);
  1148.       j:INT; loop until!(j=u.nc); -- calculate u^T.b in tmp
  1149.      s:FLT:=0.0;
  1150.      if w[j].abs>=0.000001 then    -- nonzero only if w_j is nonzero
  1151.         i:INT:=0; loop until!(i=u.nr); 
  1152.            s:=s+u[i,j]*b[i]; i:=i+1;
  1153.         end; -- loop
  1154.         s:=s/w[j];
  1155.      end; -- if
  1156.      tmp[j]:=s;
  1157.      j:=j+1;
  1158.       end; -- loop
  1159.       j:=0; loop until!(j=u.nc);
  1160.      s:FLT:=0.0;
  1161.      jj:INT:=0; loop until!(jj=u.nc);
  1162.         s:=s+v[j,jj]*tmp[jj]; jj:=jj+1
  1163.      end; -- loop
  1164.      x[j]:=s;
  1165.      j:=j+1;
  1166.       end; -- loop
  1167.    end; -- svd_back_sub
  1168.  
  1169.    inplace_linear_fit_of(vin,vout:ARRAY{VEC}):MAT
  1170.    -- Fill vin `self' to be the least squares best linear approximation
  1171.    -- relating `vin' to `vout' by: `out[i]=self.act_on(in[i])'.
  1172.    -- Return `self'.
  1173.       pre nr=vout[0].dim and nc=vin[0].dim and vout.size=vin.size is
  1174.    
  1175.       it:MAT:=create(vin.size,vin[0].dim);
  1176.       i:INT; loop until!(i=vin.size); it.inplace_row(i,vin[i]); i:=i+1; end;
  1177.       u:MAT:=create(it.nr.max(it.nc),it.nc);      
  1178.       v:MAT:=create(it.nc,it.nc); w:VEC:=w.create(it.nc);
  1179.       it.svd_in(u,w,v);
  1180.       wmax:FLT:=w.max_value; wmin:FLT:=wmax*(0.000001);
  1181.       i:=0; loop until!(i=it.nc); 
  1182.      if w[i]<=wmin then w[i]:=0.0 end; 
  1183.      i:=i+1; 
  1184.       end; -- loop
  1185.       x:VEC:=x.create(nc); b:VEC:=b.create(vout.size);
  1186.       i:=0; loop until!(i=vout[0].dim); -- get each row of self
  1187.      j:INT:=0; loop until!(j=vout.size); b[j]:=vout[j][i]; j:=j+1 end;
  1188.      it.svd_back_sub(u,w,v,b,x);
  1189.      inplace_row(i,x);
  1190.      i:=i+1;
  1191.       end; -- loop
  1192.       return self;
  1193.    end; -- inplace_linear_fit_of
  1194.  
  1195.    inplace_affine_fit_of(vin,vout:ARRAY{VEC})
  1196.    -- Fill vin `self' to be the best least squares affine map relating 
  1197.    -- `in' to `out' by: `out[i]=self.affine_act_on(vin[i])'.
  1198.       pre nr=vout[0].dim and nc=vin[0].dim+1 and vout.size=vin.size is
  1199.       it:MAT:=create(vin.size,vin[0].dim+1);
  1200.       i:INT; loop until!(i=vin.size); 
  1201.      it.inplace_row(i,vin[i]); it[i,vin[0].dim]:=1.0; -- put in affine piece
  1202.      i:=i+1; 
  1203.       end;
  1204.       u:MAT:=create(it.nr.max(it.nc),it.nc);      
  1205.       v:MAT:=create(it.nc,it.nc); w:VEC:=w.create(it.nc);
  1206.       it.svd_in(u,w,v);
  1207.       wmax:FLT:=w.max_value; wmin:FLT:=wmax*(0.000001);
  1208.       i:=0; loop until!(i=it.nc); 
  1209.      if w[i]<=wmin then w[i]:=0.0 end; 
  1210.      i:=i+1;
  1211.       end; -- loop
  1212.       x:VEC:=x.create(nc); b:VEC:=b.create(vout.size);
  1213.       i:=0; loop until!(i=vout[0].dim); -- get each row of self
  1214.      j:INT:=0; loop until!(j=vout.size); b[j]:=vout[j][i]; j:=j+1 end;
  1215.      it.svd_back_sub(u,w,v,b,x);
  1216.      inplace_row(i,x);
  1217.      i:=i+1;
  1218.       end; -- loop
  1219.    end; -- inplace_linear_fit_of
  1220.    
  1221.    inplace_weighted_linear_fit_of(vin,vout:ARRAY{VEC}, wt:ARRAY{FLT})
  1222.    -- Fill in `self' to be the least squares best linear approximation
  1223.    -- relating `vin' to `vout' by: `vout[i]=self.act_on(vin[i])'. `wt[i]'
  1224.    -- gives the weight which should be given to the ith example.
  1225.    -- (typically in `[0.,1.]' (`0.' means ignore, `1.' means full weight).
  1226.       pre
  1227.         nr=vout[0].dim and nc=vin[0].dim and vout.size=vin.size 
  1228.      and wt.size=vin.size
  1229.    is
  1230.       it:MAT:=create(vin.size,vin[0].dim);
  1231.       i:INT; loop until!(i=it.nr); it.inplace_row(i,vin[i]); i:=i+1; end;
  1232.       i:=0; loop until!(i=it.nr);    -- scale by wt
  1233.      j:INT:=0; loop until!(j=it.nc); it[i,j]:=it[i,j]*wt[i]; j:=j+1 end;
  1234.      i:=i+1;
  1235.       end; -- loop
  1236.       u:MAT:=create(it.nr.max(it.nc),it.nc);      
  1237.       v:MAT:=create(it.nc,it.nc); w:VEC:=w.create(it.nc);
  1238.       it.svd_in(u,w,v);
  1239.       wmax:FLT:=w.max_value; wmin:FLT:=wmax*(0.000001);
  1240.       i:=0; loop until!(i=it.nc); 
  1241.      if w[i]<=wmin then w[i]:=0.0 end; 
  1242.      i:=i+1;
  1243.       end; -- loop
  1244.       x:VEC:=x.create(nc); b:VEC:=b.create(vout.size);
  1245.       i:=0; loop until!(i=vout[0].dim); -- get each row of self
  1246.      j:INT:=0; loop until!(j=vout.size); b[j]:=vout[j][i]*wt[j]; j:=j+1 end;
  1247.      it.svd_back_sub(u,w,v,b,x);
  1248.      inplace_row(i,x);
  1249.      i:=i+1;
  1250.       end; -- loop
  1251.    end; -- inplace_weighted_linear_fit_of
  1252.  
  1253.    inplace_weighted_affine_fit_of(vin,vout:ARRAY{VEC}, wt:ARRAY{FLT})
  1254.    pre
  1255.         nr=vout[0].dim and nc=vin[0].dim+1 and vout.size=vin.size 
  1256.      and wt.size=vin.size
  1257.    is
  1258.       -- Fill in `self' to be the least squares best affine approximation
  1259.       -- relating `in' to `vout' by: `vout[i]=self.affine_act_on(in[i])'. 
  1260.       -- `wt[i]' gives the weight which should be given to the `i'th example. 
  1261.       -- (typically in `[0.,1.]' (`0.' means ignore, `1.' means full weight).
  1262.       it:MAT:=create(vin.size,vin[0].dim+1);
  1263.       i:INT; loop until!(i=it.nr); 
  1264.      it.inplace_row(i,vin[i]); it[i,vin[0].dim]:=wt[i];
  1265.      i:=i+1; 
  1266.       end; -- loop
  1267.       i:=0; loop until!(i=it.nr);    -- scale by wt
  1268.      j:INT:=0; loop until!(j=it.nc-1); it[i,j]:=it[i,j]*wt[i]; j:=j+1 end;
  1269.      i:=i+1;
  1270.       end; -- loop
  1271.       u:MAT:=create(it.nr.max(it.nc),it.nc);      
  1272.       v:MAT:=create(it.nc,it.nc); w:VEC:=w.create(it.nc);
  1273.       it.svd_in(u,w,v);
  1274.       wmax:FLT:=w.max_value; wmin:FLT:=wmax*(0.000001);
  1275.       i:=0; loop until!(i=it.nc); 
  1276.      if w[i]<=wmin then w[i]:=0.0 end; 
  1277.      i:=i+1; 
  1278.       end; -- loop
  1279.       x:VEC:=x.create(nc); b:VEC:=b.create(vout.size);
  1280.       i:=0; loop until!(i=vout[0].dim); -- get each row of self
  1281.      j:INT:=0; loop until!(j=vout.size); b[j]:=vout[j][i]*wt[j]; j:=j+1 end;
  1282.      it.svd_back_sub(u,w,v,b,x);
  1283.      inplace_row(i,x);
  1284.      i:=i+1;
  1285.       end; -- loop
  1286.    end; -- inplace_weighted_affine_fit_of
  1287.    
  1288. end;
  1289. -----------------------------------------------------------------------
  1290.